coef(summary(Figure.B1.1.3))["Treatment","Std. Error"], coef(summary(Figure.B1.1.4))["Treatment","Std. Error"],
coef(summary(Figure.B1.1.5))["Treatment","Std. Error"], coef(summary(Figure.B1.1.6))["Treatment","Std. Error"],
coef(summary(Figure.B1.1.7))["Treatment","Std. Error"], coef(summary(Figure.B1.1.8))["Treatment","Std. Error"],
coef(summary(Figure.B1.1.9))["Treatment","Std. Error"], coef(summary(Figure.B1.1.10))["Treatment","Std. Error"],
coef(summary(Figure.B1.1.11))["Treatment","Std. Error"], coef(summary(Figure.B1.1.12))["Treatment","Std. Error"],
coef(summary(Figure.B1.1.13))["Treatment","Std. Error"], coef(summary(Figure.B1.1.14))["Treatment","Std. Error"])
Figure.B1.1 <- data.frame(Figure.B1.Terms, Figure.B1.Estimates.1, Figure.B1.Std_Error.1)
multiplier <- 1.96
for (i in 1:nrow(Figure.B1.1)){
Figure.B1.1$ymin[i] = Figure.B1.1$Figure.B1.Estimates.1[i] - (multiplier * Figure.B1.1$Figure.B1.Std_Error.1[i])
Figure.B1.1$ymax[i] = Figure.B1.1$Figure.B1.Estimates.1[i] + (multiplier * Figure.B1.1$Figure.B1.Std_Error.1[i])
}
Figure.B1.1$Figure.B1.Terms <- factor(Figure.B1.1$Figure.B1.Terms, levels=rev(Figure.B1.1$Figure.B1.Terms), ordered=TRUE)
Figure.B1.p1 <-  ggplot(Figure.B1.1, aes(x=Figure.B1.Terms, y=Figure.B1.Estimates.1)) +
geom_hline(yintercept=0, colour="#8C2318", size=1) +  # Line at 0
geom_pointrange(aes(ymin=ymin, ymax=ymax)) +  # Ranges for each coefficient
labs(y="Coefficient of Estimates", x="", title="Support for Censorship Apparatus") +  # Labels
coord_flip() +  # Rotate the plot
theme_bw() +
theme(axis.text.x = element_text(size = 12, face = "bold"),
axis.text.y = element_text(size = 12, face = "bold")) # Nicer theme
# Heterogeneous Treatment Effect on Regime Support: Overall Satisfaction
Figure.B1.2.1 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Female == 0),])
Figure.B1.2.2 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Female == 1),])
Figure.B1.2.3 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Education %in% c(1,2)),])
Figure.B1.2.4 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Education == 3),])
Figure.B1.2.5 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Education %in% c(4,5)),])
Figure.B1.2.6 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(1,2,3)),])
Figure.B1.2.7 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(4,5)),])
Figure.B1.2.8 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(6,7,8)),])
Figure.B1.2.9 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Income %in% c(1,2)),])
Figure.B1.2.10 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Income == 3),])
Figure.B1.2.11 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Income %in% c(4,5)),])
Figure.B1.2.12 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(1,2)),])
Figure.B1.2.13 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Ideology == 3),])
Figure.B1.2.14 <- lm(Regime_Overall ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(4,5)),])
Figure.B1.2.Estimates <- c(coef(summary(Figure.B1.2.1))["Treatment","Estimate"], coef(summary(Figure.B1.2.2))["Treatment","Estimate"],
coef(summary(Figure.B1.2.3))["Treatment","Estimate"], coef(summary(Figure.B1.2.4))["Treatment","Estimate"],
coef(summary(Figure.B1.2.5))["Treatment","Estimate"], coef(summary(Figure.B1.2.6))["Treatment","Estimate"],
coef(summary(Figure.B1.2.7))["Treatment","Estimate"], coef(summary(Figure.B1.2.8))["Treatment","Estimate"],
coef(summary(Figure.B1.2.9))["Treatment","Estimate"], coef(summary(Figure.B1.2.10))["Treatment","Estimate"],
coef(summary(Figure.B1.2.11))["Treatment","Estimate"], coef(summary(Figure.B1.2.12))["Treatment","Estimate"],
coef(summary(Figure.B1.2.13))["Treatment","Estimate"], coef(summary(Figure.B1.2.14))["Treatment","Estimate"])
Figure.B1.2.Std_Error <- c(coef(summary(Figure.B1.2.1))["Treatment","Std. Error"], coef(summary(Figure.B1.2.2))["Treatment","Std. Error"],
coef(summary(Figure.B1.2.3))["Treatment","Std. Error"], coef(summary(Figure.B1.2.4))["Treatment","Std. Error"],
coef(summary(Figure.B1.2.5))["Treatment","Std. Error"], coef(summary(Figure.B1.2.6))["Treatment","Std. Error"],
coef(summary(Figure.B1.2.7))["Treatment","Std. Error"], coef(summary(Figure.B1.2.8))["Treatment","Std. Error"],
coef(summary(Figure.B1.2.9))["Treatment","Std. Error"], coef(summary(Figure.B1.2.10))["Treatment","Std. Error"],
coef(summary(Figure.B1.2.11))["Treatment","Std. Error"], coef(summary(Figure.B1.2.12))["Treatment","Std. Error"],
coef(summary(Figure.B1.2.13))["Treatment","Std. Error"], coef(summary(Figure.B1.2.14))["Treatment","Std. Error"])
Figure.B1.2 <- data.frame(Figure.B1.Terms, Figure.B1.2.Estimates, Figure.B1.2.Std_Error)
for (i in 1:nrow(Figure.B1.2)){
Figure.B1.2$ymin[i] = Figure.B1.2$Figure.B1.2.Estimates[i] - (multiplier * Figure.B1.2$Figure.B1.2.Std_Error[i])
Figure.B1.2$ymax[i] = Figure.B1.2$Figure.B1.2.Estimates[i] + (multiplier * Figure.B1.2$Figure.B1.2.Std_Error[i])
}
Figure.B1.2$Figure.B1.Terms <- factor(Figure.B1.2$Figure.B1.Terms, levels=rev(Figure.B1.2$Figure.B1.Terms), ordered=TRUE)
Figure.B1.p2 <- ggplot(Figure.B1.2, aes(x=Figure.B1.Terms, y=Figure.B1.2.Estimates)) +
geom_hline(yintercept=0, colour="#8C2318", size=1) +  # Line at 0
geom_pointrange(aes(ymin=ymin, ymax=ymax)) +  # Ranges for each coefficient
labs(y="Coefficient of Estimates", x="", title="Regime Support: Overall") +  # Labels
coord_flip() +  # Rotate the plot
theme_bw() +
theme(axis.text.x = element_text(size = 12, face = "bold"),
axis.text.y = element_text(size = 12, face = "bold")) # Nicer theme
# Heterogeneous Treatment Effect on Regime Support: Central Government
Figure.B1.3.1 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Female == 0),])
Figure.B1.3.2 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Female == 1),])
Figure.B1.3.3 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Education %in% c(1,2)),])
Figure.B1.3.4 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Education == 3),])
Figure.B1.3.5 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Education %in% c(4,5)),])
Figure.B1.3.6 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(1,2,3)),])
Figure.B1.3.7 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(4,5)),])
Figure.B1.3.8 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(6,7,8)),])
Figure.B1.3.9 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Income %in% c(1,2)),])
Figure.B1.3.10 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Income == 3),])
Figure.B1.3.11 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Income %in% c(4,5)),])
Figure.B1.3.12 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(1,2)),])
Figure.B1.3.13 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Ideology == 3),])
Figure.B1.3.14 <- lm(Regime_Central ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(4,5)),])
Figure.B1.3.Estimates <- c(coef(summary(Figure.B1.3.1))["Treatment","Estimate"], coef(summary(Figure.B1.3.2))["Treatment","Estimate"],
coef(summary(Figure.B1.3.3))["Treatment","Estimate"], coef(summary(Figure.B1.3.4))["Treatment","Estimate"],
coef(summary(Figure.B1.3.5))["Treatment","Estimate"], coef(summary(Figure.B1.3.6))["Treatment","Estimate"],
coef(summary(Figure.B1.3.7))["Treatment","Estimate"], coef(summary(Figure.B1.3.8))["Treatment","Estimate"],
coef(summary(Figure.B1.3.9))["Treatment","Estimate"], coef(summary(Figure.B1.3.10))["Treatment","Estimate"],
coef(summary(Figure.B1.3.11))["Treatment","Estimate"], coef(summary(Figure.B1.3.12))["Treatment","Estimate"],
coef(summary(Figure.B1.3.13))["Treatment","Estimate"], coef(summary(Figure.B1.3.14))["Treatment","Estimate"])
Figure.B1.3.Std_Error <- c(coef(summary(Figure.B1.3.1))["Treatment","Std. Error"], coef(summary(Figure.B1.3.2))["Treatment","Std. Error"],
coef(summary(Figure.B1.3.3))["Treatment","Std. Error"], coef(summary(Figure.B1.3.4))["Treatment","Std. Error"],
coef(summary(Figure.B1.3.5))["Treatment","Std. Error"], coef(summary(Figure.B1.3.6))["Treatment","Std. Error"],
coef(summary(Figure.B1.3.7))["Treatment","Std. Error"], coef(summary(Figure.B1.3.8))["Treatment","Std. Error"],
coef(summary(Figure.B1.3.9))["Treatment","Std. Error"], coef(summary(Figure.B1.3.10))["Treatment","Std. Error"],
coef(summary(Figure.B1.3.11))["Treatment","Std. Error"], coef(summary(Figure.B1.3.12))["Treatment","Std. Error"],
coef(summary(Figure.B1.3.13))["Treatment","Std. Error"], coef(summary(Figure.B1.3.14))["Treatment","Std. Error"])
Figure.B1.3 <- data.frame(Figure.B1.Terms, Figure.B1.3.Estimates, Figure.B1.3.Std_Error)
for (i in 1:nrow(Figure.B1.3)){
Figure.B1.3$ymin[i] = Figure.B1.3$Figure.B1.3.Estimates[i] - (multiplier * Figure.B1.3$Figure.B1.3.Std_Error[i])
Figure.B1.3$ymax[i] = Figure.B1.3$Figure.B1.3.Estimates[i] + (multiplier * Figure.B1.3$Figure.B1.3.Std_Error[i])
}
Figure.B1.3$Figure.B1.Terms <- factor(Figure.B1.3$Figure.B1.Terms, levels=rev(Figure.B1.3$Figure.B1.Terms), ordered=TRUE)
Figure.B1.p3 <- ggplot(Figure.B1.3, aes(x=Figure.B1.Terms, y=Figure.B1.3.Estimates)) +
geom_hline(yintercept=0, colour="#8C2318", size=1) +  # Line at 0
geom_pointrange(aes(ymin=ymin, ymax=ymax)) +  # Ranges for each coefficient
labs(y="Coefficient of Estimates", x="", title="Regime Support: Central") +  # Labels
coord_flip() +  # Rotate the plot
theme_bw() +
theme(axis.text.x = element_text(size = 12, face = "bold"),
axis.text.y = element_text(size = 12, face = "bold")) # Nicer theme
# Heterogeneous Treatment Effect on Regime Support: Local Government
Figure.B1.4.1 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Female == 0),])
Figure.B1.4.2 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Female == 1),])
Figure.B1.4.3 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Education %in% c(1,2)),])
Figure.B1.4.4 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Education == 3),])
Figure.B1.4.5 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Education %in% c(4,5)),])
Figure.B1.4.6 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(1,2,3)),])
Figure.B1.4.7 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(4,5)),])
Figure.B1.4.8 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(6,7,8)),])
Figure.B1.4.9 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Income %in% c(1,2)),])
Figure.B1.4.10 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Income == 3),])
Figure.B1.4.11 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Income %in% c(4,5)),])
Figure.B1.4.12 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(1,2)),])
Figure.B1.4.13 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Ideology == 3),])
Figure.B1.4.14 <- lm(Regime_Local ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(4,5)),])
Figure.B1.4.Estimates <- c(coef(summary(Figure.B1.4.1))["Treatment","Estimate"], coef(summary(Figure.B1.4.2))["Treatment","Estimate"],
coef(summary(Figure.B1.4.3))["Treatment","Estimate"], coef(summary(Figure.B1.4.4))["Treatment","Estimate"],
coef(summary(Figure.B1.4.5))["Treatment","Estimate"], coef(summary(Figure.B1.4.6))["Treatment","Estimate"],
coef(summary(Figure.B1.4.7))["Treatment","Estimate"], coef(summary(Figure.B1.4.8))["Treatment","Estimate"],
coef(summary(Figure.B1.4.9))["Treatment","Estimate"], coef(summary(Figure.B1.4.10))["Treatment","Estimate"],
coef(summary(Figure.B1.4.11))["Treatment","Estimate"], coef(summary(Figure.B1.4.12))["Treatment","Estimate"],
coef(summary(Figure.B1.4.13))["Treatment","Estimate"], coef(summary(Figure.B1.4.14))["Treatment","Estimate"])
Figure.B1.4.Std_Error <- c(coef(summary(Figure.B1.4.1))["Treatment","Std. Error"], coef(summary(Figure.B1.4.2))["Treatment","Std. Error"],
coef(summary(Figure.B1.4.3))["Treatment","Std. Error"], coef(summary(Figure.B1.4.4))["Treatment","Std. Error"],
coef(summary(Figure.B1.4.5))["Treatment","Std. Error"], coef(summary(Figure.B1.4.6))["Treatment","Std. Error"],
coef(summary(Figure.B1.4.7))["Treatment","Std. Error"], coef(summary(Figure.B1.4.8))["Treatment","Std. Error"],
coef(summary(Figure.B1.4.9))["Treatment","Std. Error"], coef(summary(Figure.B1.4.10))["Treatment","Std. Error"],
coef(summary(Figure.B1.4.11))["Treatment","Std. Error"], coef(summary(Figure.B1.4.12))["Treatment","Std. Error"],
coef(summary(Figure.B1.4.13))["Treatment","Std. Error"], coef(summary(Figure.B1.4.14))["Treatment","Std. Error"])
Figure.B1.4 <- data.frame(Figure.B1.Terms, Figure.B1.4.Estimates, Figure.B1.4.Std_Error)
for (i in 1:nrow(Figure.B1.4)){
Figure.B1.4$ymin[i] = Figure.B1.4$Figure.B1.4.Estimates[i] - (multiplier * Figure.B1.4$Figure.B1.4.Std_Error[i])
Figure.B1.4$ymax[i] = Figure.B1.4$Figure.B1.4.Estimates[i] + (multiplier * Figure.B1.4$Figure.B1.4.Std_Error[i])
}
Figure.B1.4$Figure.B1.Terms <- factor(Figure.B1.4$Figure.B1.Terms, levels=rev(Figure.B1.4$Figure.B1.Terms), ordered=TRUE)
Figure.B1.p4 <- ggplot(Figure.B1.4, aes(x=Figure.B1.Terms, y=Figure.B1.4.Estimates)) +
geom_hline(yintercept=0, colour="#8C2318", size=1) +  # Line at 0
geom_pointrange(aes(ymin=ymin, ymax=ymax)) +  # Ranges for each coefficient
labs(y="Coefficient of Estimates", x="", title="Regime Support: Local") +  # Labels
coord_flip() +  # Rotate the plot
theme_bw() +
theme(axis.text.x = element_text(size = 12, face = "bold"),
axis.text.y = element_text(size = 12, face = "bold")) # Nicer theme
# Heterogeneous Treatment Effect on Willingness to Protest
Figure.B1.5.1 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Female == 0),])
Figure.B1.5.2 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Female == 1),])
Figure.B1.5.3 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Education %in% c(1,2)),])
Figure.B1.5.4 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Education == 3),])
Figure.B1.5.5 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Education %in% c(4,5)),])
Figure.B1.5.6 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(1,2,3)),])
Figure.B1.5.7 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(4,5)),])
Figure.B1.5.8 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Age_Group %in% c(6,7,8)),])
Figure.B1.5.9 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Income %in% c(1,2)),])
Figure.B1.5.10 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Income == 3),])
Figure.B1.5.11 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Income %in% c(4,5)),])
Figure.B1.5.12 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(1,2)),])
Figure.B1.5.13 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Ideology == 3),])
Figure.B1.5.14 <- lm(Regime_Protest ~ Treatment, data = Experiment[which(Experiment$Ideology %in% c(4,5)),])
Figure.B1.5.Estimates <- c(coef(summary(Figure.B1.5.1))["Treatment","Estimate"], coef(summary(Figure.B1.5.2))["Treatment","Estimate"],
coef(summary(Figure.B1.5.3))["Treatment","Estimate"], coef(summary(Figure.B1.5.4))["Treatment","Estimate"],
coef(summary(Figure.B1.5.5))["Treatment","Estimate"], coef(summary(Figure.B1.5.6))["Treatment","Estimate"],
coef(summary(Figure.B1.5.7))["Treatment","Estimate"], coef(summary(Figure.B1.5.8))["Treatment","Estimate"],
coef(summary(Figure.B1.5.9))["Treatment","Estimate"], coef(summary(Figure.B1.5.10))["Treatment","Estimate"],
coef(summary(Figure.B1.5.11))["Treatment","Estimate"], coef(summary(Figure.B1.5.12))["Treatment","Estimate"],
coef(summary(Figure.B1.5.13))["Treatment","Estimate"], coef(summary(Figure.B1.5.14))["Treatment","Estimate"])
Figure.B1.5.Std_Error <- c(coef(summary(Figure.B1.5.1))["Treatment","Std. Error"], coef(summary(Figure.B1.5.2))["Treatment","Std. Error"],
coef(summary(Figure.B1.5.3))["Treatment","Std. Error"], coef(summary(Figure.B1.5.4))["Treatment","Std. Error"],
coef(summary(Figure.B1.5.5))["Treatment","Std. Error"], coef(summary(Figure.B1.5.6))["Treatment","Std. Error"],
coef(summary(Figure.B1.5.7))["Treatment","Std. Error"], coef(summary(Figure.B1.5.8))["Treatment","Std. Error"],
coef(summary(Figure.B1.5.9))["Treatment","Std. Error"], coef(summary(Figure.B1.5.10))["Treatment","Std. Error"],
coef(summary(Figure.B1.5.11))["Treatment","Std. Error"], coef(summary(Figure.B1.5.12))["Treatment","Std. Error"],
coef(summary(Figure.B1.5.13))["Treatment","Std. Error"], coef(summary(Figure.B1.5.14))["Treatment","Std. Error"])
Figure.B1.5 <- data.frame(Figure.B1.Terms, Figure.B1.5.Estimates, Figure.B1.5.Std_Error)
for (i in 1:nrow(Figure.B1.5)){
Figure.B1.5$ymin[i] = Figure.B1.5$Figure.B1.5.Estimates[i] - (multiplier * Figure.B1.5$Figure.B1.5.Std_Error[i])
Figure.B1.5$ymax[i] = Figure.B1.5$Figure.B1.5.Estimates[i] + (multiplier * Figure.B1.5$Figure.B1.5.Std_Error[i])
}
Figure.B1.5$Figure.B1.Terms <- factor(Figure.B1.5$Figure.B1.Terms, levels=rev(Figure.B1.5$Figure.B1.Terms), ordered=TRUE)
Figure.B1.p5 <- ggplot(Figure.B1.5, aes(x=Figure.B1.Terms, y=Figure.B1.5.Estimates)) +
geom_hline(yintercept=0, colour="#8C2318", size=1) +  # Line at 0
geom_pointrange(aes(ymin=ymin, ymax=ymax)) +  # Ranges for each coefficient
labs(y="Coefficient of Estimates", x="", title="Willingness to Protest") +  # Labels
coord_flip() +  # Rotate the plot
theme_bw() +
theme(axis.text.x = element_text(size = 12, face = "bold"),
axis.text.y = element_text(size = 12, face = "bold")) # Nicer theme
Figure.B1 <- grid.arrange(Figure.B1.p1,
Figure.B1.p2,
Figure.B1.p3,
Figure.B1.p4,
Figure.B1.p5, nrow = 3)
# --------------------------------------------------
# Table B6
# --------------------------------------------------
# Outcome variables
formula <- c("Censor_Support ~ ", "Regime_Overall ~ ", "Regime_Central ~ ",
"Regime_Local ~ ", "Regime_Protest ~ ")
# Study 1 Results (Difference in Means)
lm_list_s1 <- lapply(formula, function(f) lm(paste0(f, "Treatment"), data = Study1))
names(lm_list_s1) <- c("lm_censor_s1", "lm_overall_s1", "lm_central_s1", "lm_local_s1", "lm_protest_s1")
# Study 2 Results (Difference in Means)
lm_list_s2 <- lapply(formula, function(f) lm(paste0(f, "Treatment"), data = Study2))
names(lm_list_s2) <- c("lm_censor_s2", "lm_overall_s2", "lm_central_s2", "lm_local_s2", "lm_protest_s2")
# Two Studies Combined Results (Difference in Means)
lm_list_c <- lapply(formula, function(f) lm(paste0(f, "Treatment"), data = Experiment))
names(lm_list_c) <- c("lm_censor_c", "lm_overall_c", "lm_central_c", "lm_local_c", "lm_protest_c")
# Combining Models into a list
models <- c(lm_list_s1, lm_list_s2, lm_list_c)
# Point Estimates
Estimates <- round(sapply(models, function(x) coef(summary(x))["Treatment","Estimate"]),3)
# p-values
p.value <- paste("[", round(sapply(models, function(x) coef(summary(x))["Treatment","Pr(>|t|)"]),5), "]")
# Adjusting p-values using Benjamini-Hochberg p-value correction
Adj.p.value <- paste("[", round(p.adjust(sapply(models, function(x) coef(summary(x))["Treatment","Pr(>|t|)"]), method = "BH"),5), "]")
table.B6 <- rbind(Estimates[c(1,6,11)], p.value[c(1,6,11)], Adj.p.value[c(1,6,11)],
Estimates[c(2,7,12)], p.value[c(2,7,12)], Adj.p.value[c(2,7,12)],
Estimates[c(3,8,13)], p.value[c(3,8,13)], Adj.p.value[c(3,8,13)],
Estimates[c(4,9,14)], p.value[c(4,9,14)], Adj.p.value[c(4,9,14)],
Estimates[c(5,10,15)], p.value[c(5,10,15)], Adj.p.value[c(5,10,15)])
colnames(table.B6) <- c("Study 1", "Study 2", "Combined")
rownames(table.B6) <- c("Support for Censorship Apparatus", "   p-value", "   adjusted p-value",
"Overall Satisfaction of China", "   p-value", "   adjusted p-value",
"Assessment of Central Governmen", "   p-value", "   adjusted p-value",
"Assessment of Local Government", "   p-value", "   adjusted p-value",
"Willingness to Protest", "   p-value", "   adjusted p-value")
table.B6
# --------------------------------------------------
# Figure B2
# --------------------------------------------------
# Create a new variable for Explicit Support for Censorship Apparatus (Binary)
Study2 <- Study2 %>%
mutate(across(
c(Censor_Support),
list(Binary = ~ case_when(
. %in% c(4,5) ~ 1,
. %in% c(1,2,3) ~ 0,
TRUE ~ NA_real_
))
))
## Group Means for Explicit Censorship Support
Figure.B2 <-  Study2 %>%
drop_na(Treatment) %>%
group_by(Treatment) %>%
get_summary_stats(Censor_Support_Binary, type = "mean_sd") %>%
mutate(variable = case_when(
variable == "Censor_Support_Binary" ~ "Explicit Support"
),
Treatment = case_when(
Treatment == 1 ~ "Treatment",
Treatment == 0 ~ "Control"
),
# Standard errors for group mean of Explicit Censorship Support
se = sd / sqrt(n),
# Confidence interval for group mean Explicit Censorship Support
ymin = mean - 1.96 * se,
ymax = mean + 1.96 * se)
# Calculating Implicit Censorship Support
# Implicit support in the treatment group
lmt <- lm(Censor_Support_Implicit ~ Censor_Support_LE_T, data = Study2[Study2$Treatment == 1,])
# Implicit support in the control group
lmc <- lm(Censor_Support_Implicit ~ Censor_Support_LE_T, data = Study2[Study2$Treatment == 0,])
# Create a data frame for Figure B2
# Data frame for implicit support
Treatment <- c("Control", "Treatment")
Variable <- c(rep("Implicit Support",2))
Estimates <- c(coef(summary(lmc))["Censor_Support_LE_T","Estimate"],
coef(summary(lmt))["Censor_Support_LE_T","Estimate"])
Std_Error <- c(coef(summary(lmc))["Censor_Support_LE_T","Std. Error"],
coef(summary(lmt))["Censor_Support_LE_T","Std. Error"])
lower <- c(confint(lmc, level = 0.95)["Censor_Support_LE_T",][1],
confint(lmt, level = 0.95)["Censor_Support_LE_T",][1])
upper <- c(confint(lmc, level = 0.95)["Censor_Support_LE_T",][2],
confint(lmt, level = 0.95)["Censor_Support_LE_T",][2])
Figure.B2.2  <- data.frame(Treatment, Variable, Estimates, Std_Error, lower, upper)
# Combining the data frame for explicit support and the data frame for implicit support
Figure.B2  <- Figure.B2[,c(1,2,4,6,7,8)]
colnames(Figure.B2) <- colnames(Figure.B2.2)
Figure.B2 <- rbind(Figure.B2, Figure.B2.2)
# Plot out implicit and explicit support
pd = position_dodge2(width = 0.3, reverse = TRUE)
ggplot(Figure.B2, aes(x=Treatment, y = Estimates, shape = Variable)) +
# Point estimates
geom_point(position = pd, size = 3) +
# Confidence intervals
geom_linerange(aes(x = Treatment,
ymin = lower,
ymax = upper),
size = 0.5,
position = pd, show.legend = FALSE) +
ggtitle("Support for Censorship Apparatus") +
labs(y="Proportion of Respondents", x="") +  # Labels
theme_bw() +
scale_y_continuous(limits = c(0.2,0.8)) +
theme(axis.text = element_text(size = 12, face = "bold"),
axis.title = element_text(size=12),
axis.text.y = element_text(size = 12),
legend.title = element_blank()) # Legend label for shape
rm(list=ls())
library("caret")
library("irr")
library("data.table")
# Load the data
TableD1.data <- fread("TableD1.csv")
TableE1.data <- fread("TableE1.csv")
TableE2.data <- fread("TableE2.csv")
TableE3.data <- fread("TableE3.csv")
# Create a confusion matrix
confusion_matrix <- confusionMatrix(data = factor(TableD1.data$Coder2, levels = unique(TableD1.data$Coder2)),
reference = factor(TableD1.data$Coder1, levels = unique(TableD1.data$Coder1)))
# Calculate precision, recall, and F1-score by category
Table.D1 <- round(confusion_matrix$byClass[, c("Precision", "Recall", "F1")], 2)
# Calculate macro precision, recall, and F1-score
Macro <- round(c(mean(Table.D1[,"Precision"]),
mean(Table.D1[,"Recall"]),
mean(Table.D1[,"F1"])),2)
# Sort the rows alphabetically by category name
Table.D1 <- as.data.frame(t(Table.D1[order(rownames(Table.D1)), ]))
# Modify the column names
colnames(Table.D1) <- substr(colnames(Table.D1),
nchar(colnames(Table.D1)) - 2,
nchar(colnames(Table.D1)))
# Attach Macro precision, recall, and F1-score
Table.D1$Macro <- Macro
Table.D1
coder.confusion <- table(TableD1.data$Coder1, TableD1.data$Coder2)
coder.accuracy <- paste(round(sum(diag(coder.confusion))/sum(coder.confusion),4)*100, "%", sep = "")
coder.accuracy
TableD1.data$Coder1_Political <- ifelse(TableD1.data$Coder1 %in% c("COL", "CRI", "GOV"),1,0)
TableD1.data$Coder2_Political <- ifelse(TableD1.data$Coder2 %in% c("COL", "CRI", "GOV"),1,0)
coder.confusion.pol <- table(TableD1.data$Coder1_Political, TableD1.data$Coder2_Political)
coder.accuracy.pol <- paste(round(sum(diag(coder.confusion.pol))/sum(coder.confusion.pol),4)*100, "%", sep = "")
coder.accuracy.pol
round(kappa2(TableD1.data)$value,2)
# Create a function for calculating macro F1 score
f1_score <- function(actual, predicted) {
confusion_matrix <- table(actual, predicted)
precisions <- diag(confusion_matrix) / rowSums(confusion_matrix)
recalls <- diag(confusion_matrix) / colSums(confusion_matrix)
f1_scores <- 2 * (precisions * recalls) / (precisions + recalls)
macro_f1 <- mean(f1_scores, na.rm = TRUE)
return(macro_f1)
}
round(kappa2(TableD1.data)$value,2)
View(Table.D1)
View(TableD1.data)
round(kappa2(TableD1.data[,1:2])$value,2)
# Create a function for calculating macro F1 score
f1_score <- function(actual, predicted) {
confusion_matrix <- table(actual, predicted)
precisions <- diag(confusion_matrix) / rowSums(confusion_matrix)
recalls <- diag(confusion_matrix) / colSums(confusion_matrix)
f1_scores <- 2 * (precisions * recalls) / (precisions + recalls)
macro_f1 <- mean(f1_scores, na.rm = TRUE)
return(macro_f1)
}
# Calculate the macro F1 score for each model prediction
Table.E1 <- data.frame(Macro_F1_Score = round(c(
f1_score(TableE1.data$Original_Label, TableE1.data$BERT_Prediction),
f1_score(TableE1.data$Original_Label, TableE1.data$LogitRegression_Prediction),
f1_score(TableE1.data$Original_Label, TableE1.data$PaLM_Prediction),
f1_score(TableE1.data$Original_Label, TableE1.data$XGBoost_Prediction),
f1_score(TableE1.data$Original_Label, TableE1.data$RandomForest_Prediction),
f1_score(TableE1.data$Original_Label, TableE1.data$VotingClassifier_Prediction),
f1_score(TableE1.data$Original_Label, TableE1.data$DecisionTree_Prediction),
f1_score(TableE1.data$Original_Label, TableE1.data$NeuralNetwork_Prediction),
f1_score(TableE1.data$Original_Label, TableE1.data$Word2Vec_Prediction)),4))
rownames(Table.E1) <- c("Fine-tuned Pre-Train Chinese BERT",
"Logistic Regression (Ridge)",
"Pattern Learning and Matching (PaLM)",
"Extreme Gradient Boosting (XGBoost)",
"Random Forest",
"Ensemble Classifier (Voting)",
"Decision Tree",
"Neural Network",
"Word2Vec Embedding")
Table.E1
# Create a confusion matrix
confusion_matrix.E2 <- confusionMatrix(data = factor(TableE2.data$Outsample_Prediction, levels = unique(TableE2.data$Outsample_Prediction)),
reference = factor(TableE2.data$Coder1, levels = unique(TableE2.data$Coder1)))
# Calculate precision, recall, and F1-score by category
Table.E2 <- round(confusion_matrix.E2$byClass[, c("Precision", "Recall", "F1")], 2)
# Calculate macro precision, recall, and F1-score
Macro.E2 <- round(c(mean(Table.E2[,"Precision"]),
mean(Table.E2[,"Recall"]),
mean(Table.E2[,"F1"])),2)
# Sort the rows alphabetically by category name
Table.E2 <- as.data.frame(t(Table.E2[order(rownames(Table.E2)), ]))
# Modify the column names
colnames(Table.E2) <- substr(colnames(Table.E2),
nchar(colnames(Table.E2)) - 2,
nchar(colnames(Table.E2)))
# Attach Macro precision, recall, and F1-score
Table.E2$Macro <- Macro.E2
Table.E2
Insample.confusion <- table(TableE2.data$Coder1, TableE2.data$Insample_Prediction)
Insample.accuracy <- round(sum(diag(Insample.confusion))/sum(Insample.confusion),2)
Insample.accuracy
weights <- 1/(sum(1/table(TableE2.data$Coder1))) * 1/table(TableE2.data$Coder1)
weighted_precision <- round(sum(as.numeric(Table.E2["Precision",1:9]) * weights),2)
weighted_precision
weighted_recall  <- round(sum(as.numeric(Table.E2["Recall",1:9]) * weights),2)
weighted_recall
weighted_macro_F1 <- (weighted_recall + weighted_precision)/2
weighted_macro_F1
TableE2.data$Coder_Political <- ifelse(TableE2.data$Coder1 %in% c("COL", "CRI", "GOV"),1,0)
TableE2.data$Outsample_Political <- ifelse(TableE2.data$Outsample_Prediction %in% c("COL", "CRI", "GOV"),1,0)
confusion_pol <- table(TableE2.data$Outsample_Political, TableE2.data$Coder_Political)
weights_pol <- 1/(sum(1/table(TableE2.data$Coder_Political))) * 1/table(TableE2.data$Coder_Political)
weighted_precision_pol <- round(sum(diag(confusion_pol) / rowSums(confusion_pol) * weights_pol),2)
weighted_precision_pol
weighted_recall_pol  <- round(sum(diag(confusion_pol) / colSums(confusion_pol) * weights_pol),2)
weighted_recall_pol
weighted_macro_F1_pol <- round((weighted_recall_pol + weighted_precision_pol)/2,2)
weighted_macro_F1_pol
# Calculate the proportion of each specific category
TableE3.proportion <- round(table(TableE3.data$Multinomial_Prediction)/nrow(TableE3.data),4)*100
Table.E3 <- data.frame(
General_Category = c(rep("Highly Political",4),rep("Moderately Political",3),rep("Non-Political",5)),
Specific_Category = c("Collective Action", "Govt Criticism", "Other Govt-related", "Total",
"Business", "Foreign", "Total",
"Entertainment", "Advertisement", "Culture", "Others", "Total"),
Logistical_Regression_Ridge = c(TableE3.proportion["COL"],TableE3.proportion["CRI"],TableE3.proportion["GOV"],
sum(TableE3.proportion[c("COL", "CRI", "GOV")]),
TableE3.proportion["BET"],TableE3.proportion["FOR"],
sum(TableE3.proportion[c("BET", "FOR")]),
TableE3.proportion["ESX"],TableE3.proportion["ADS"],TableE3.proportion["LCT"],TableE3.proportion["OTH"],
sum(TableE3.proportion[c("ESX", "ADS", "LCT", "OTH")]))
)
Table.E3[,"Logistical_Regression_Ridge"] <- paste(Table.E3[,"Logistical_Regression_Ridge"], "%", sep = "")
Table.E3
rm(list=ls())
# Packages ----
library("dplyr")
library("ggplot2")
library("data.table")
# Load the data
FreeWeChat <- fread("FreeWeChatResults.csv")
WeChatScope <- fread("WeChatScopeResults.csv")
WeiboScope <- fread("WeiboScopeResults.csv")
# Raw count of each category divided by total number of observations
# FreeWeChat censorship data
FreeWeChat_Table <- data.frame(
Laebl = colnames(FreeWeChat[,3:11]),
Category =  c("Advertisement", "Business", "Collective Action",
"Govt Criticism", "Entertainment", "Foreign",
"Other Govt-related", "Culture", "Others"),
General_Category = c("Non-Political", "Moderately-Political", "Highly-Political",
"Highly-Political", "Non-Political", "Moderately-Political",
"Highly-Political", "Non-Political", "Non-Political"),
Sum = colSums(FreeWeChat)[3:11],
Percentage = round((colSums(FreeWeChat)[3:11] / sum(colSums(FreeWeChat)[3:11])) * 100,2)
)
# WeChatScope censorship data
WeChatScope_Table <- data.frame(
Laebl = colnames(WeChatScope[,3:11]),
Category =  c("Advertisement", "Business", "Collective Action",
"Govt Criticism", "Entertainment", "Foreign",
"Other Govt-related", "Culture", "Others"),
General_Category = c("Non-Political", "Moderately-Political", "Highly-Political",
"Highly-Political", "Non-Political", "Moderately-Political",
"Highly-Political", "Non-Political", "Non-Political"),
Sum = colSums(WeChatScope)[3:11],
Percentage = round((colSums(WeChatScope)[3:11] / sum(colSums(WeChatScope)[3:11])) * 100,2)
)
# WeiboScope censorship data
WeiboScope_Table <- data.frame(
Laebl = colnames(WeiboScope[,3:11]),
Category =  c("Advertisement", "Business", "Collective Action",
"Govt Criticism", "Entertainment", "Foreign",
"Other Govt-related", "Culture", "Others"),
General_Category = c("Non-Political", "Moderately-Political", "Highly-Political",
"Highly-Political", "Non-Political", "Moderately-Political",
"Highly-Political", "Non-Political", "Non-Political"),
Sum = colSums(WeiboScope)[3:11],
Percentage = round((colSums(WeiboScope)[3:11] / sum(colSums(WeiboScope)[3:11])) * 100,2)
)
# Merging the three censorship data source
Merged_Table <- merge(merge(WeChatScope_Table[,c(2,3,5)],
FreeWeChat_Table[,c(2,3,5)],
by = c("Category","General_Category"), all = TRUE),
WeiboScope_Table[,c(2,3,5)],
by = c("Category","General_Category"), all = TRUE)
colnames(Merged_Table)[3:5] <- c("WeChatScope", "FreeWeChat", "WeiboScope")
Merged_Table <- Merged_Table[,c(2,1,3:5)]
Table.1 <- rbind(as.character(Merged_Table[Merged_Table$Category == "Collective Action",]),
as.character(Merged_Table[Merged_Table$Category == "Govt Criticism",]),
as.character(Merged_Table[Merged_Table$Category == "Other Govt-related",]),
c("Highly-Political", "Total",
colSums(Merged_Table[Merged_Table$General_Category == "Highly-Political",3:5])),
as.character(Merged_Table[Merged_Table$Category == "Business",]),
as.character(Merged_Table[Merged_Table$Category == "Foreign",]),
c("Moderately-Political", "Total",
colSums(Merged_Table[Merged_Table$General_Category == "Moderately-Political",3:5])),
as.character(Merged_Table[Merged_Table$Category == "Entertainment",]),
as.character(Merged_Table[Merged_Table$Category == "Advertisement",]),
as.character(Merged_Table[Merged_Table$Category == "Culture",]),
as.character(Merged_Table[Merged_Table$Category == "Others",]),
c("Non-Political", "Total",
colSums(Merged_Table[Merged_Table$General_Category == "Non-Political",3:5])))
Table.1[,3:5] <- paste(Table.1[,3:5], "%", sep = "")
colnames(Table.1) <- colnames(Merged_Table)
Table.1
# Combine highly political categories and combine other categories
## Calculate the proportion of political vs. non-political categories
FreeWeChat$Political <- rowSums(FreeWeChat[, c("COL", "GOV", "CRI")]) / rowSums(FreeWeChat[, 3:11])
